data = readtable('../../data/Netherlands/Holland_PF_v2.xlsx');
date = data.Var1;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T = length(date(startdate:enddate));

option.detrend = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

time = data.Var1(startdate:enddate);
pi = data.pi(startdate:enddate);
rpcgdpgr = data.x(startdate:enddate);
tau = data.tau(startdate:enddate);
g = data.g(startdate:enddate);
ynom10 = data.y10(startdate:enddate);
surplus = tau - g;

debt_book = data.debt_gdp(startdate:enddate); % book value of debt/gdp ratio
debt = data.debt_gdp_mkt(startdate:enddate);

%%convenience yield
cy = 0.015;
ynom10 = ynom10 + cy;
ynom10 = log(ynom10 + 1);

tau = tau + debt * cy;

deltalogg = [log(data.g(startdate)) - log(data.g(startdate - 1)); log(g(2:end)) - log(g(1:end - 1))];
deltalogtau = [log(data.tau(startdate)) - log(data.tau(startdate - 1)); log(tau(2:end)) - log(tau(1:end - 1))];

x0 = mean(rpcgdpgr);
pi0 = mean(pi);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
y0nom_10 = mean(ynom10);
surplus0 = mean(surplus);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if option.detrend == 1
    tts(1) = tau(1);
    gts(1) = g(1);

    for t = 2:T
        gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
        tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
    end

    disp('mean surplus')
    100 * [mean(tts - gts) mean(surplus)]
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define the ordering in the VAR
inflpos = 1;
ypos = 2;
gdppos = 3;
dtpos = 4;
coint_tax = 5;
dgpos = 6;
coint_spending = 7;

%% raw data
clear Xraw;
Xraw(:, inflpos) = pi;
Xraw(:, gdppos) = rpcgdpgr;
Xraw(:, ypos) = ynom10;
Xraw(:, dtpos) = deltalogtau;
Xraw(:, coint_tax) = log(tau);
Xraw(:, dgpos) = deltalogg;
Xraw(:, coint_spending) = log(g);

% VAR elements
infl = pi - pi0;
x = rpcgdpgr - x0;
yield10 = ynom10 - y0nom_10;
dg = deltalogg - g0;
dt = deltalogtau - tau0;

X2(:, inflpos) = infl;
X2(:, ypos) = yield10;
X2(:, gdppos) = x;
X2(:, dtpos) = dt;
X2(:, coint_tax) = log(t) - mean(log(t));
X2(:, dgpos) = dg;
X2(:, coint_spending) = log(g) - mean(log(g));

% trend adjustment
if option.detrend == 1
    X2(:, coint_tax) = log(tts) - mean(log(tts));
    X2(:, coint_spending) = log(gts) - mean(log(gts));
end

z_var = X2(2:end, :);
z_varlag = X2(1:end - 1, :);
Y = z_var;
X = z_varlag;

N = cols(X2);
I = eye(N);
I_pi = I(:, inflpos);
I_gdp = I(:, gdppos);
I_y10 = I(:, ypos);
I_dt = I(:, dtpos);
I_dg = I(:, dgpos);
I_coint_tax = I(:, coint_tax);
I_coint_spending = I(:, coint_spending);

% actual, non-detrended series of tax and spending
X2t = X2;
X2t(:, coint_tax) = log(tau) - mean(log(tau));
X2t(:, coint_spending) = log(g) - mean(log(g));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Psi and Sig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi = zeros(N, N);
Psi_se = zeros(N, N);

for i = [inflpos ypos gdppos dtpos dgpos]
    regr = ols(Y(:, i), [ones(T - 1, 1), X(:, :)]);
    c(i) = regr.beta(1);
    c_se(i) = regr.bstd(1);
    Psi(i, :) = regr.beta(2:end)';
    Psi_se(i, :) = regr.bstd(2:end)';
    R2(i) = regr.rsqr;
    R2bar(i) = regr.rbar;
    eps(:, i) = Y(:, i) - c(i) - X * Psi(i, :)';

    clear regr;
end

Psi(coint_tax, :) = Psi(dtpos, :);
Psi(coint_spending, :) = Psi(dgpos, :);

Psi(coint_tax, coint_tax) = Psi(coint_tax, coint_tax) + 1;
Psi(coint_spending, coint_spending) = Psi(coint_spending, coint_spending) + 1;

Psi_se(coint_tax, :) = Psi_se(dtpos, :);
Psi_se(coint_spending, :) = Psi_se(dgpos, :);

Sigma = cov(eps(:, [inflpos ypos gdppos dtpos dgpos]));
tmp = chol(Sigma, 'lower');
Sig = zeros(N, N);
Sig([inflpos ypos gdppos], [inflpos ypos gdppos]) = ...
    tmp(1:3, 1:3);

Sig(dtpos, [inflpos ypos gdppos dtpos]) = tmp(4, 1:4);
Sig(dgpos, [inflpos ypos gdppos dtpos dgpos]) = tmp(5, 1:5);

Sig(coint_tax, :) = Sig(dtpos, :);
Sig(coint_spending, :) = Sig(dgpos, :);

eps = [eps(:, [inflpos ypos gdppos]), ...
                eps(:, dtpos), eps(:, dtpos), ...
                eps(:, dgpos), eps(:, dgpos)];
Sigma = Sig * Sig';

tstat = Psi ./ Psi_se;
